#include "Element.h"

int main(int argc, char* argv[]){
  /*
  Element<double,2> E;
  E.set_n_dofs(3);
  E.set_accuracy(2);
  Point<double,2> P1({2.3,4.5},23,0);
  Point<double,2> P2({0.0,-1.9},24,1);
  E.set_global_dofs(0,P1);
  E.set_global_dofs(1,P2);
  std::cout << "n_dofs: " << E.get_n_dofs() << "," << std::endl;
  std::cout << "global: " << E.get_global_dofs(0).get_coord() << ", " <<E.get_global_dofs(1).get_coord()
	    <<  std::endl;
  std::cout << "global index: " << E.get_global_dofs(0).get_index() << ", " <<E.get_global_dofs(1).get_index()
	    <<  std::endl;
  std::cout << "global boundary: " << E.get_global_dofs(0).get_boundary_mark() << ", " <<E.get_global_dofs(1).get_boundary_mark()
	    <<  std::endl;
  std::cout << "accracy: " << E.get_accuracy() << "," << std::endl;
  */
  TD_P1Element<double> Ele;
  /*
  std::cout << "n_dofs: " << Ele.get_n_dofs() << "," << std::endl;
  std::cout << "local coord: " << Ele.get_local_dofs(0).get_coord() << ", " <<Ele.get_local_dofs(1).get_coord()
    << ", " << Ele.get_local_dofs(2).get_coord() <<  std::endl;
  */

  Point<double,2> local_P({0.5,0.5},-1,1);
  Point<double,2> global_t1({1.0,0},12,1);
  Point<double,2> global_t2({2.0,0},13,1);
  Point<double,2> global_t3({1.5,2},14,0);
  Ele.set_global_dofs(0,global_t1);
  Ele.set_global_dofs(1,global_t2);
  Ele.set_global_dofs(2,global_t3);
  /*
  Point<double,2> local_result = Ele.local_to_global(local_P);
  std::cout << "result: " << std::endl;
  std::cout << "coord: " << local_result.get_coord() << std::endl;
  std::cout << "index: " << local_result.get_index() << std::endl;
  */
  /*
  Point<double,2> P1({0.2,0.3},-1,1);
  double X0 = Ele.local_basis_function(0,P1);
  double X1 = Ele.local_basis_function(1,P1);
  double X2 = Ele.local_basis_function(2,P1);
  std::cout << "local basis function: " << X0 << "," << X1 << "," << X2 << std::endl;
  Vector<double,2> G0 = Ele.local_basis_gradient(0);
  Vector<double,2> G1 = Ele.local_basis_gradient(1);
  Vector<double,2> G2 = Ele.local_basis_gradient(2);
  std::cout << "local basis gradient: " << G0 << "," << G1 << "," << G2 << std::endl;
  */
  /*
  Matrix<double,2,2> M = Ele.global_to_local_Jacobi();
  std::cout << "g2l Jacobi:" << std::endl << M << std::endl;
  double det = Ele.global_to_local_det();
  std::cout << "g2l Jacobi det:" << det << std::endl;
  Matrix<double,2,2> rM = Ele.local_to_global_Jacobi();
  std::cout << "l2g Jacobi:" << std::endl << rM << std::endl;
  */
  int N = 3;
  Ele.read_quad_info("data/triangle.tmp_geo",N);
  std::cout << "acc: " << Ele.get_accuracy() << std::endl;
  const std::valarray<Point<double,2> >& Quad_P = Ele.get_quad_points();
  const std::valarray<double>& Quad_w = Ele.get_quad_weights();
  int n_quad = Quad_P.size();
  /*
  for (int i = 0 ; i < n_quad ; i++){
    std::cout << i << " : coord=" << std::endl << Quad_P[i].get_coord() << std::endl
	      << "index=" << Quad_P[i].get_index() << std::endl;
    std::cout << "weight=" << Quad_w[i] << std::endl;
  }
  */
  for (int i = 0 ; i < n_quad ; i++){
    std::cout << i << " : coord=" << std::endl << Ele.get_quad_points(i).get_coord() << std::endl
	      << "index=" << Ele.get_quad_points(i).get_index() << std::endl;
    std::cout << "weight=" << Ele.get_quad_weights(i) << std::endl;
  }
  return 0;
}
